Cite this Article
Shan Ming-Lei, Zhu Chang-Ping, Yao Cheng, Yin Cheng, Jiang Xiao-Yan. Pseudopotential multi-relaxation-time lattice Boltzmann model for cavitation bubble collapse with high density ratio. Chinese Physics B, 2016, 25(10): 104701  
Permissions
Pseudopotential multi-relaxation-time lattice Boltzmann model for cavitation bubble collapse with high density ratio
Shan Ming-Lei1, 2, Zhu Chang-Ping1, 2, †, , Yao Cheng1, Yin Cheng1, Jiang Xiao-Yan1
College of Internet of Things Engineering, Hohai University, Changzhou 213022, China
Jiangsu Key Laboratory of Power Transmission and Distribution Equipment Technology, Hohai University, Changzhou 213022, China

 

† Corresponding author. E-mail: cpzhu5126081@163.com

Project supported by the National Natural Science Foundation of China (Grant Nos. 11274092 and 1140040119) and the Natural Science Foundation of Jiangsu Province, China (Grant No. SBK2014043338).

Abstract
Abstract

The dynamics of the cavitation bubble collapse is a fundamental issue for the bubble collapse application and prevention. In the present work, the modified forcing scheme for the pseudopotential multi-relaxation-time lattice Boltzmann model developed by Li Q et al. [Li Q, Luo K H and Li X J 2013 Phys. Rev. E 87 053301] is adopted to develop a cavitation bubble collapse model. In the respects of coexistence curves and Laplace law verification, the improved pseudopotential multi-relaxation-time lattice Boltzmann model is investigated. It is found that the thermodynamic consistency and surface tension are independent of kinematic viscosity. By homogeneous and heterogeneous cavitation simulation, the ability of the present model to describe the cavitation bubble development as well as the cavitation inception is verified. The bubble collapse between two parallel walls is simulated. The dynamic process of a collapsing bubble is consistent with the results from experiments and simulations by other numerical methods. It is demonstrated that the present pseudopotential multi-relaxation-time lattice Boltzmann model is applicable and efficient, and the lattice Boltzmann method is an alternative tool for collapsing bubble modeling.

PACS: 47.11.Qr;47.55.Ca;47.55.dd;47.55.dp
Keyword:lattice Boltzmann method;pseudopotential model;bubble collapse;improved forcing scheme
1. Introduction

Cavitation is a unique phenomenon in liquid, and has made extensive research efforts in science and engineering. Due to the various dynamic effects during the cavitation generation and cavitation bubble collapse, such as bubble wall oscillation and the formations of shockwave and micro-jet, the potential applications of cavitation have been explored and are being extensively studied in the fields of medicine, biology, lab on chip, sewage treatment, and material surface cleaning treatment.[1–5] However, as too many phenomena are involved, a theoretical model of cavitation bubble collapse is difficult to be established, and for particular boundary conditions, the analytical solution is even impossible. Therefore, the numerical simulation becomes a powerful way to gain a better understanding. The conventional numerical simulation methods of cavitation bubble mainly include the finite volume method (FVM),[6–8] the finite element method (FEM),[9] and the boundary element method (BEM).[10] In these partial differential equation based numerical simulations, the methods to track or capture the interfaces (such as the volume of fluid (VOF)[11] or level set method (LSM)[12]) are required, which are often demanding computationally. In addition, the Poisson equation needs to be solved to satisfy the continuity equation, which drastically reduces the computational efficiency.[13]

In recent years, the lattice Boltzmann method (LBM), which is a mesoscopic approach based on the kinetic Boltzmann equation, has emerged as a powerful tool for simulating multiphase flow problems.[14–18] As an alternative tool for the numerical simulations and investigations of multiphase flows, the LBM provides many advantages including clear physical pictures, easy implementation of boundary conditions, and fully parallel algorithms.[19] Generally, the existing LBM multiphase models can be mostly classified as being in four categories: the color-gradient method,[20] the pseudopotential method (or Shan-Chen model),[21,22] the free-energy method,[23,24] and the phase-field method.[25] The pseudopotential method is widely and successfully used in the LBM multiphase community due to its conceptual simplicity and computation efficiency. In the pseudopotential method, the fluid interactions are mimicked by an interparticle potential, from which a non-monotomic equation of state (EOS) can be obtained. As a result, the separation of fluid phases or components can be achieved automatically in this method, and the methods to track or capture the interfaces are not required yet. Moreover, the pressure can be calculated from EOS efficiently instead of the Poisson equation.

Recently, the pseudopotential LBM was firstly introduced into the issue of cavitation by Sukop and Or.[26] In the following several years, research efforts have been made to investigate the mechanism of cavitation. Chen et al.[27] simulated the cavitating bubble growth by using the modified pseudopotential LB model with the EDM force scheme. The results in quiescent flows agree fairly well with the solution of the Rayleigh–Plesset equation. Mishra et al.[28] introduced a model of cavitation based on the pseudopotential LB model that allows for coupling between the hydrodynamics of a collapsing cavity and supported solute chemical species. Using a pseudopotential LB model, Daemi et al.[29] modeled the bubble cluster in an acoustic field, and verified that the deformation and coalescence phenomena in bubble cluster dynamics can be simulated by LBM. However, the above-mentioned research on cavitation did not include the cavitation bubble collapse process, which is essential for the study of the cavitation phenomenon, and the density ratio between liquid phase and gas phase was limited below 100.

For the pseudopotential LB model, the density ratio is a synthetic problem relating to the lattice model, forcing scheme, thermodynamic consistency, numerical stability etc., and great efforts have been devoted to this issue.[30–35] In the present work, the modified forcing scheme for pseudopotential Multi-Relaxation-Time (MRT) LB developed by Li Q et al.[33] is adopted to develop a cavitation bubble collapse model. Essentially, Li et al.’s scheme is an approximate approach to thermodynamic consistency by adjusting forcing scheme. This approach retains the concise and high efficiency of the original LBM, and achieves the stable LB multiphase model with a large density ratio.

The rest of the present paper is organized as follows. Section 2 will briefly introduce the pseudopotential MRT-LB model and the forcing scheme modified by Li et al. The numerical analysis of the modified pseudopotential MRT-LB model for a bubble will be given in Section 3. In Section 4, the validity and feasibility of the method for cavitation and bubble collapse will be verified. Finally, a brief conclusion will be made in Section 5.

2. Modified pseudopotential MRT-LB model
2.1. Pseudopotential MRT-LB model

The pseudopotential LB model, also known as the Shan–Chen model, was developed by Shan and Chen in 1993.[21] In the pseudopotential method, the fluid interactions are mimicked by an interparticle potential, which is now widely called peseudopotential. In the original pseudopotential LB model, the single-relaxation-time (SRT) collision operator was employed. In recent years, the MRT collision operator has been verified, showing that it is superior to the SRT operator in terms of numerical stability. The MRT-LB evolution equation can be given as follows:[33]

where fα is the density distribution function, is its equilibrium distribution, t is the time, x is the spatial position, eα is the discrete velocity along the α-th direction, δt is the time step, is the forcing term in the velocity space, M is an orthogonal transformation matrix, Λ a diagonal matrix, and for D2Q9 lattice, Λ is given by

Through the transformation matrix M, f, and feq can be projected onto the moment space via m = M f and meq = M feq, and the collision step of the MRT-LB equation (1) can be rewritten as[33]

where I is the unit tensor, and S is the forcing term in the moment space with (I − 0.5Λ)S = M F′. For D2Q9 lattice, meq can be given by

where is the macroscopic density, v is the macroscopic velocity and . Then the streaming step of the MRT-LB equation can be formulated as

where f* = M−1m*. The macroscopic velocity in Eq. (4) is calculated by

where F = (Fx,Fy) for two-dimensional (2D) space is the force acting on the system.

For the pseudopotential LB model, the F in Eq. (6) is given by[36]

where ψ (x) is the interparticle potential, G is the interaction strength, and w1−4 = 1/3 and w5−8 = 1/12 are the weights, F is incorporated via S with a specific forcing scheme. For the MRT-LB method, the usual forcing scheme for D2Q9 lattice can be given by

2.2. Li et al.’s forcing scheme for pseudopotential MRT-LB model

For the pseudopotential multiphase LB model, the vapor density ρg, liquid density ρl and the pressure p0 in an equilibrium coexistence fluid, should satisfy the mechanical stability condition expressed as[33,36]

where ψ′ = dψ/ρ, ɛ is a parameter related with the interaction range.[36] pEOS−PL in Eq. (9) is the pressure from the non-ideal EOS of the pseudopotential LB model, which is given by

where cs is the lattice sound speed. Using Eqs. (9) and (10), the coexistence curve of the pseudopotential LB model can be obtained by numerical integration. From the thermodynamic theory view, the equilibrium coexistence fluid should also satisfy the Maxwell equal-area rule which determines the thermodynamic coexistence and is expressed as:

Here, pEOS−TD is the EOS in the thermodynamic theory. Without special treatment, the mechanical stability condition determined by Eq. (9) will deviate from the solution given by the Maxwell construction of Eq. (10). That is usually called the thermodynamic inconsistency in the pseudopotential LB model.

To address the thermodynamic inconsistency, the integrands in Eqs. (9) and (11) should be equal or approximate to each other. Each integrand is naturally divided into two parts, i.e., the bracketed part which relates directly to EOS and the outside part which relates directly to ρ. So, two different strategies have been developed to solve the thermodynamic inconsistency. The first strategy is to guarantee that the outside parts in Eqs. (9) and (11) are equal to each other,[37] i.e., ψ′/ψ1+ɛ = 1/ρ2. In this strategy framework, the specific form of ψ is proposed as follows:[37]

However, it is difficult for the first strategy to make pEOS−PL = pEOS−TD. In addition, the form EOS given by Eq. (10) is fixed. So, the adoption of different EOS is impossible.

The second strategy is to guarantee that the bracketed parts in Eqs. (9) and (10) are equal to each other at first[38] (i.e., pEOS−PL = pEOS−TD = pEOS), then to reduce the deviation between ψ′/ψ1+ɛ and 1/ρ2. As a result, ψ can be given by

Here, G is used to ensure that the whole term inside the square root is positive.[38] The second strategy allows different EOS to be involved, and is widely used in pseudopotential LB models.[15]

In the second strategy framework, Li et al.[33] proposed an MRT version forcing scheme to achieve thermodynamic consistency. For the D2Q9 lattice, Li et al.’s forcing scheme can be given by

Compared with Eq. (8), S1 and S2, the components of S, are modified with a term which can be tuned by ɛ. Through theoretical analysis, Li et al. found that there exists a suitable ɛ which can make the mechanical stability solution approximately identical to the solution given by the thermodynamic consistency requirement in a wide range of temperature.[32,33] Li et al.’s forcing scheme has the following advantages: (i) maintaining a uniform layout with a general form of the LB forcing scheme; (ii) achieving thermodynamic consistency only by tuning one constant parameter; and (c) fully retaining the LBM’s advantages of being simple and efficient.

3. Numerical simulations and analyses

In this section, numerical investigations for the bubble will be conducted with Li et al.’s improved pseudopotential LB model. Firstly, the improved model will be investigated for bubble simulations from three aspects: the optimal ɛ for bubble simulations, the effect on coexistence curves from relaxation time τυ, and the validation of the Laplace law. Then, the homogeneous and heterogeneous cavitation will be simulated to investigate the feasibility of the presented LB model for the cavitation phenomenon.

3.1. Investigation of improved model

Three numerical simulations of bubbles are considered to investigate Li et al.’s improved pseudopotential MRT-LB model. The first one is to simulate the stationary bubble, which can be used to obtain a numerical coexistence curve approximate to that given by the Maxwell construction by tuning ɛ. The second numerical simulation is to investigate the effect of relaxation time. The last simulation is conducted to validate the law of Laplace. In the present work, the Carnahan–Starling (CS) EOS is adopted, which can be given by[38]

where and b = 0.18727RTc/pc with Tc and pc being the critical temperature and pressure, respectively. In our simulations, we set a = 0.5, b = 1, R = 1, and δt = 1. A 201 × 201 lattice is adopted in the simulations of this section. A bubble with radius r0 = 70 is initially placed in the center of the domain. The density field is initialized as[31]

where (x0,y0) is the center of the domain, W is the prescribed width of the phase interface and it is set to be 5 in the present work, ‘tanh’ is the hyperbolic tangent function and tanh (x) = (e2x−1)/(e2x+1). The periodical boundary conditions are applied to the x and y directions. Due to the use of CS EOS, ψ adopts the form in Eq. (13), and G = −1 is used. The relaxation times are chosen as follows: , , and .

Li et al. have estimated that the optimal value of ɛ should be between 1 and 2.[32] For the stationary liquid droplet with radius r0 = 50, they considered that ɛ = 1.76 is an optimal value. With this value, the coexistence curves obtained from LBM simulations are in good agreement with those given by the Maxwell equal-area construction. For a stationary vapor bubble, the coexistence curves of the case τυ = 0.6 with different ɛ values are shown in Fig. 1. From the figure we can see that when ɛ = 1.0 the vapor branch of the coexistence curve significantly deviates from the solution of the Maxwell equal-area construction and the achievable lowest reduced temperature is around 0.7Tc. For ɛ = 2.0, the LB model works well at 0.5Tc, but there are still obvious deviations for the coexistence curves. In contrast, the coexistence curve at ɛ = 1.86 agrees quite well with the solution of the Maxwell equal-area construction in a wide range of temperature. The maximum density ratio is beyond 700, which is close to the density ratio between water and vapor in nature. Figure 1 demonstrates that for a vapor bubble, the improved MRT pseduopotential LB model is capable of achieving thermodynamic consistency, a large density ratio and stability by tuning ɛ. In addition, it should be noted that the optimal value of ɛ for the bubble case is larger than that for the droplet case. This may be associated with the effect of surface tension and the difference in compressibility between vapor phase and liquid phase.

Fig. 1. Comparison of the numerical coexistence curves by LBM with different ɛ values (τυ = 0.6).

Furthermore, the effect of relaxation time τυ, which is associated with kinematic viscosity, is investigated. The coexistence curves of the cases τυ = 0.51, τυ = 0.6, and τυ = 0.8 as ɛ = 1.86 are shown in Fig. 2. We can see that their good agreement can also be achieved. It demonstrates that thermodynamic consistency is independent of kinematic viscosity for Li et al.’s improved MRT pseduopotential LB model.

Fig. 2. Numerical coexistence curves by LBM with different τυ values (ɛ = 1.86).

Satisfying the Laplace law is an important benchmark test for multiphase flows. The surface tension force, meanwhile, can be calculated when the Laplace law is verified. For the bubble case, the Laplace law can be given as

where pout and pin are the fluid pressures outside and inside the bubble, respectively, σ is the surface tension, and r is the radius of the bubble. The results at different reduced temperatures are shown in Fig. 3 for the cases τυ = 0.6 and τυ = 0.8. From Fig. 3, the linear relationship between Δp and 1/r can be clearly observed. The values of line slope (σ) are almost the same for the cases of τυ = 0.6 and τυ = 0.8, and decrease with temperature increasing. It indicates that the improved MRT pseduopotential LB model satisfies the Laplace law, and the independence between σ and τυ makes it more convenient to investigate the physical mechanism of the multiphase flows.

Fig. 3. Numerical coexistence curves by LBM with different τυ values (ɛ = 1.86).
3.2. Cavitation inception and development

The inception of a cavitation bubble is essentially a phase transition from liquid to vapor. It happens when the pure liquid is exposed to such a high negative pressure that makes the density of liquid in the range of coexistence density.

In the present work, the cavitation inception phenomenon is investigated by the proposed MRT pseduopotential LB model. A (lx,ly) = (201,201) lattice is adopted. The whole density field is initialized as ρinit = 0.2609 which is between the critical density ρc = 0.1305 and the liquid density ρl = 0.4541 at 0.5Tc with τυ = 0.6 and ɛ = 1.86. Then a density perturbation of ρinit/108 is set at (ly − 1)/2. The periodical boundary conditions are applied to x and y directions. The densities at ((lx − 1)/2,ly − 1) and ((lx − 1)/2, (ly − 1)/2) are detected. The density evolutions with time are shown in Fig. 4. In our simulation, the cavitation occurs as a linear band because the initial density perturbation is linear and the boundary conditions are periodical. At the first stage, the liquid phase and the vapor phase are not separated apparently. At about 180 time steps, the vapor phase is firstly separated at (ly − 1)/2. After severe density fluctuations, the density at ((lx − 1)/2, (ly − 1)/2) is stable near the density of vapor phase ρg = 6.269× 10−4. In contrast, the density at ((lx − 1)/2, ly − 1) is fluctuated for much more time before approaching to the density of the liquid phase. Finally, a stable fluid system including liquid phase and vapor phase is formed. These simulation results demonstrate that the proposed MRT pseduopotential LB model is applicable to describing the process of cavitation inception essentially.

Fig. 4. Density evolution in the process of cavitation inception, simulated by LBM.

If there is a cavitation nucleus in the static fluid, the later development of the cavitation bubble is determined by the critical radius of cavitation nucleus, which can be given as follows:[39]

where r* is the critical radius of the cavitation bubble, σ is the surface tension and Δp is the pressure difference between vapor and liquid phases relative to the flat interface. σ can be obtained from the Laplace law verification in Subsection 3.1, and the surface tension at 0.5Tc is used here, i.e., σ = 0.01051. Specific Δp can be achieved by pressure boundaries[26] or by setting the density of liquid domain artificially. The second approach is used here. By tuning the liquid phase, a critical radius of 11.3 is achieved finally. In this simulation, a 401×401 lattice is adopted and the periodical boundary conditions are applied to x and y directions. For verifying the effect of cavitation nuclei critical radius, the cavitation bubbles with radii greater and smaller than critical radius are initialized in the center of the simulation domain, respectively. The evolutions of the bubbles are shown in Fig. 5. We can see that a bubble with radius just below the critical value cannot overcome the energy barrier and eventually condenses, while a bubble with a radius just above the critical value expands gradually and reaches a new equilibrium eventually. It demonstrates that the proposed MRT pseduopotential LB model is valid to simulate the cavitation bubble development as well as the cavitation inception.

Fig. 5. Effect of the cavitation nucleus critical radius, verified by LBM.
4. Bubble collapse between two parallel solid walls

A set of investigations for the proposed MRT pseduopotential LB model demonstrates that the LB model can establish cavitation numerical models in line with the law of physics. And with the stability including such a high density ratio close to that between water and vapor in reality, the LB model is promising to establish an efficient cavitation numerical model. In this section, we will establish a numerical model to investigate the mechanism of the bubble collapse between two parallel solid walls.

The dynamics of cavitation bubble near a solid wall is an important issue for understanding the mechanism of surface damages in fluid machinery. It has been widely studied and most of these studies are based on the semi-infinite fluid domain.[40–42] However, in the case of the bubble between two parallel solid walls, the collapse of the bubble will be significantly different from in the semi-infinite assumption.[3,43] Although several investigations have been conducted based on macroscopic methods, the LB modeling of cavitation bubble collapse between two solid walls is still in its infancy. In this section, using the proposed MRT pseduopotential LB model with high density ratio, we make an attempt to simulate bubble collapse between two parallel solid walls in the 2D case.

4.1. Computational domain

The computational domain for the bubble collapse between two parallel solid walls is shown in Fig. 6. In the present simulation, a 1801×301 lattice is adopted. Using Eq. (16), a spherical bubble with radius R0 is initialized in a stationary fluid between two parallel walls. H is the distance between upper and lower walls, i.e., H = 301. And d is the distance between the bubble center and the lower wall. Pv and P∞ represent the pressure of vapor inside the bubble and the ambient pressure outside the bubble, respectively.

Fig. 6. Schematic diagram of the computational domain.

For describing the initial simulation states quantitatively, two dimensionless parameters are introduced as follows:

where λ is defined as the stand-off distance between the bubble center and the lower wall; κ is defined as the relative location of the bubble to the channel between two walls. In this simulation, κ is set to be 0.5, i.e., the bubble center is located on the horizontal centerline of the channel; λ is set to be 0.92 according to Ref. [44]. So the parameters of computational domain can be assigned, i.e., d = 150 and R0 = 138.

In order to simulate the bubble collapse process, a positive pressure difference ΔP = P∞ − Pv is achieved by artificially tuning the initial liquid density based on the equilibrium state. The periodical boundary condition is applied to the x direction, and a no-slip bounce-back boundary condition is implemented in the y direction. The parameters of the LB model are set to be T/Tc = 0.5, τυ = 0.51 and ɛ = 1.86.

4.2. Simulation results

The simulation results of bubble collapse between two parallel walls by the proposed MRT pseduopotential LB model are shown in Fig. 7, where the experimental results[44] and the VOF simulations[43] are also presented for comparison. In the whole collapse process, the longitudinal length of the bubble is almost fixed due to the retarding effect of two parallel walls. In contrast, the bubble laterally shrinks step by step, presenting a dumbbell type, and then splits into two bubbles after the collapse. Finally, the jets towards upper and lower walls are formed. This dynamic process is consistent with the experimental result and VOF simulation result.

Fig. 7. Comparison of results among (a) experiments, (b) LBM simulations, (c) VOF simulations.

To investigate the dynamic characteristics in more details when the bubble collapses and the jet forms, the pressure distribution, the regional density and velocity distributions at these two key moments are shown in Figs. 8 and 9. The dumbbell-type bubble splits into two sub-bubbles at the moment of collapse. From Fig. 8, we can find that the impetus of collapse comes from the horizontal liquid flowing to the bubble neck. A colliding interface can be found obviously at the tips of the two cone-shaped sub-bubbles. As a result, a higher pressure/density region forms at the bubble neck. This higher pressure/density region becomes the leading source power forming the jets.

Fig. 8. Pressure distribution (left) and the regional density and velocity distributions (right) as bubble collapses.
Fig. 9. Pressure distribution (left) and the regional density and velocity distributions (right) as a jet forms.

Due to the regional higher pressure at the tips of sub-bubbles, the tips are flatted and then sag. Accordingly, there emerges the high flow velocity towards the insides of two sub-bubbles, respectively. Then the jets of the collapsing bubble form, as shown in Fig. 9. The high velocity jet is considered to be the leading factor in wall damage. The detailed investigations of the interaction mechanisms between the high velocity jet of the collapsing bubble and the parallel walls will be developed by the present LBM model in our future work.

5. Conclusions

For the modeling of a collapsing bubble, an improved MRT pseduopotential LB model is investigated in some respects such as thermodynamic consistency, Laplace law, homogeneous and heterogeneous cavitation. Then the bubble collapse between two parallel walls is simulated. The dynamic process of the collapsing bubble is consistent with the results from experiments and simulations obtained by other numerical methods.

The improved forcing scheme developed by Li et al. provides a convenient and efficient approach to achieving thermodynamic consistency. Adopting Li et al.’s forcing scheme, we find that the thermodynamic consistency and surface tension are independent of kinematic viscosity, which makes it superior to investigate the physical mechanism of multiphase flows. Even at 0.5Tc with τυ = 0.51, the improved MRT pseduopotential LB model also has enough stability to describe the collapsing bubble with a high density ratio beyond 700. It is demonstrated that the present MRT pseduopotential LB model is applicable and efficient, and the LBM is an alternative tool for collapsing bubble modeling.

Reference
1Brennen C E2005Fundamentals of Multiphase FlowNew YorkCambridge University Press192219–22
2Hashmi AYu GReilly-Collette MHeiman GXu J 2012 Lab Chip 12 4216
3Azam F IKarri BOhl S WKlaseboer EKhoo B C 2013 Phys. Rev. 88 043006
4Xu ZYasuda KLiu X J 2015 Chin. Phys. 24 104301
5van Wijngaarden L 2016 Ultrason. Sonochem. 29 524
6Chau S WHsu K LKouh J SChen Y J 2004 J. Mar. Sci. Technol. 8 147
7Cai JLi XLiu BHuai X 2014 Chin. Sci. Bull. 59 1580
8Liu BCai JHuai X LLi F C 2013 J. Heat Trans-T. ASME 136 022901
9Hsu C YLiang C CNguyen A TTeng T L 2014 Ocean. Eng. 81 29
10Zhang A MLiu Y L 2015 J. Comput. Phys. 294 208
11Samiei EShams MEbrahimi R 2011 Eur. J. Mech. B-Fluid 30 41
12Sussman M 2003 J. Comput. Phys. 187 110
13He Y LWang YLi Q2009Lattice Boltzmann method: theory and applicationsBeijingScience Press1101–10(in Chinese)
14Succi S2001The lattice Boltzmann equation for fluid dynamics and beyondOxfordClarendon405040–50
15Huang HSukop MLu X2015Multiphase Lattice Boltzmann Methods: Theory and ApplicationOxfordJohn Wiley & Sons131–3
16Li QLuo K HKang Q JHe Y LChen QLiu Q 2016 Prog. Energ. Combust. 52 62
17Xu A GZhang G CLi Y JLi H2014Prog. Phys.34136(in Chinese)
18Chen LKang QMu YHe Y LTao W Q 2014 Int. J. Heat Mass. Trans. 76 210
19Chen S YDoolen G D 1998 Ann. Rev. Fluid Mech. 30 329
20Gunstensen A KRothman D HZaleski SZanetti G 1991 Phys. Rev. 43 4320
21Shan X WChen H D 1993 Phys. Rev. 47 1815
22Shan XChen H 1994 Phys. Rev. 49 2941
23Swift M ROsborn W RYeomans J M 1995 Phys. Rev. Lett. 75 830
24Swift M ROrlandini EOsborn WYeomans J 1996 Phys. Rev. 54 5041
25He X YChen S YZhang R Y 1999 J. Comput. Phys. 152 642
26Sukop MOr D 2005 Phys. Rev. 71 046703
27Chen X PZhong C WYuan X L 2011 Comput. Math. Appl. 61 3577
28Mishra S KDeymier P AMuralidharan KFrantziskonis GPannala SSimunovic S 2010 Ultrason. Sonochem. 17 258
29Daemi MTaeibi-Rahni MMassah H 2015 Chin. Phys. 24 024302
30Kupershtokh A LMedvedev D AKarpov D I 2009 Comput. Math. Appl. 58 965
31Huang H BKrafczyk MLu X 2011 Phys. Rev. 84 046710
32Li QLuo K HLi X J 2012 Phys. Rev. 86 016709
33Li QLuo K HLi X J 2013 Phys. Rev. 87 053301
34Xie H QZeng ZZhang L Q 2016 Chin. Phys. 25 014702
35Song B WRen FHu H BHuang Q G 2015 Chin. Phys. 24 014703
36Shan X W 2008 Phys. Rev. 77 066702
37Sbragaglia MShan X W 2011 Phys. Rev. 84 036703
38Yuan PSchaefer L 2006 Phys. Fluids 18 042101
39Or DTuller M2002Water Resour. Res.3819
40Plesset M SChapman R B 1971 J. Fluid. Mech. 47 283
41Brujan E AKeen G SVogel ABlake J R 2002 Phys. Fluids 14 85
42Lauer EHu X YHickel SAdams N A 2012 Comput. Fluids 69 1
43Liu BCai JHuai X L 2014 Int. J. Heat Mass. Trans. 78 830
44Ishida HNuntadusit CKimoto HNakagawa TYamamoto T2001CAV 2001: Fourth International Symposium on CavitationJune 20–23, 2001Pasadena, USAA5 003